﻿using System;
using System.Collections.Generic;
using System.Drawing;
using System.Drawing.Imaging;
using System.Runtime.InteropServices;
using System.Windows.Forms;
using MentalAlchemy.Atomics;
using MentalAlchemy.Molecules.Image.Struct;

namespace MentalAlchemy.Molecules.Image
{
	/// <summary>
	/// Class to unite classes and methods useful for image analysis.
	/// Includes calculation of different image parameters such as
	/// - gradients,
	/// - decomposition to color channels,
	/// - various stats (min/max/mean brightness, gray levels frequencies etc.),
	/// - salient points (using Harris and Harris-Laplace detectors.),
	/// - image quality (using Munteanu-Rosa and modified Munteanu-Rosa criteria).
	/// </summary>
	public class ImageAnalysis
	{
		#region - Gradient calculation. -
		public static float[,] CalculateGradient(Bitmap image, Orientation orient)
		{
			int width = image.Width, height = image.Height;
			var res = new float[height, width];

			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			var stride = bmpData.Stride;
			var alignByte = stride - width * 3;
			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			if (orient == Orientation.Horizontal)
			{	// calculate gradient along X-axis.
				var count = 0;
				alignByte += 3;
				for (int i = 0; i < height; ++i)
				{
					count += 3;	// skip the first byte.
					for (int j = 1; j < width - 1; ++j, count += 3)
					{
						res[i, j] = bytes[count + 3] - bytes[count - 3];
					}
					count += alignByte;
				}

				for (int i = 0; i < height; ++i)
				{
					res[i, width - 1] = 0.0f;
				}
			}
			else
			{	// calculate gradient along Y-axis.
				int count = 0;
				count += stride;	// skip the first line.
				for (int i = 1; i < height - 1; ++i)
				{
					for (int j = 0; j < width; ++j, count += 3)
					{
						res[i, j] = bytes[count + stride] - bytes[count - stride];
					}
					count += alignByte;
				}

				for (int i = 0; i < width; ++i)
				{
					res[height - 1, i] = 0.0f;
				}
			}

			image.UnlockBits(bmpData);

			return res;
		}
		#endregion

		#region - Color channels. -
		/// <summary>
		/// 
		/// </summary>
		/// <param name="image"></param>
		/// <param name="colorBitmaps"></param>
		public static void GetColorChannels(Bitmap image, ref Bitmap[] colorBitmaps)
		{
			if (image != null && (image.PixelFormat == PixelFormat.Format24bppRgb
								|| image.PixelFormat == PixelFormat.Format32bppArgb
								|| image.PixelFormat == PixelFormat.Format32bppRgb))
			{
				colorBitmaps = new Bitmap[3];
				colorBitmaps[0] = new Bitmap(image.Size.Width, image.Size.Height, PixelFormat.Format24bppRgb);
				colorBitmaps[1] = new Bitmap(image.Size.Width, image.Size.Height, PixelFormat.Format24bppRgb);
				colorBitmaps[2] = new Bitmap(image.Size.Width, image.Size.Height, PixelFormat.Format24bppRgb);

				var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
				var bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

				int bytesSize = image.Height * bmpData.Stride;
				var bytes = new byte[bytesSize];
				var bytesR = new byte[bytesSize];
				var bytesG = new byte[bytesSize];
				var bytesB = new byte[bytesSize];
				Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

				int width = image.Width * 3,
					alignByte = bmpData.Stride - width,
					count = 0;
				for (var i = 0; i < bytesSize; )
				{
					var limit = i + width;
					for (var j = i; j < limit; j += 3, count += 3)
					{
						bytesB[count] = bytesB[count + 1] = bytesB[count + 2] = bytes[j];
						bytesG[count] = bytesG[count + 1] = bytesG[count + 2] = bytes[j + 1];
						bytesR[count] = bytesR[count + 1] = bytesR[count + 2] = bytes[j + 2];
					}
					i = limit + alignByte;
					count += alignByte;
				}

				image.UnlockBits(bmpData);

				// set R layer
				var bmpDataR = colorBitmaps[0].LockBits(rect, ImageLockMode.WriteOnly, colorBitmaps[0].PixelFormat);
				Marshal.Copy(bytesR, 0, bmpDataR.Scan0, bytesSize);
				colorBitmaps[0].UnlockBits(bmpDataR);

				// set G layer
				var bmpDataG = colorBitmaps[1].LockBits(rect, ImageLockMode.WriteOnly, colorBitmaps[1].PixelFormat);
				Marshal.Copy(bytesG, 0, bmpDataG.Scan0, bytesSize);
				colorBitmaps[1].UnlockBits(bmpDataG);

				// set B layer
				var bmpDataB = colorBitmaps[2].LockBits(rect, ImageLockMode.WriteOnly, colorBitmaps[2].PixelFormat);
				Marshal.Copy(bytesB, 0, bmpDataB.Scan0, bytesSize);
				colorBitmaps[2].UnlockBits(bmpDataB);
			}
		}
		#endregion

		#region - Image stats calculation. -

		#region - Global stats. -
		/// <summary>
		/// [molecule]
		/// 
		/// Calculates min brightness for the given image.
		///  Note! The input image is supposed to be a grayscale one!
		/// </summary>
		/// <param name="image">Input image.</param>
		/// <returns>Min brightness value</returns>
		public static float CalculateMinBrightness(Bitmap image)
		{
			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);

			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			var minB = 0;
			for (var i = 0; i < bytes.Length; i += 3)
			{
				if (minB > bytes[i])
				{
					minB = bytes[i];
				}
			}
			image.UnlockBits(bmpData);

			return minB;
		}

		/// <summary>
		/// [molecule]
		/// 
		/// Calculates max brightness for the given image.
		/// Note! The input image is supposed to be a grayscale one!
		/// </summary>
		/// <param name="image">Input image.</param>
		/// <returns>Max brightness value</returns>
		public static float CalculateMaxBrightness(Bitmap image)
		{
			//float res = 0.0;

			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);

			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			var maxB = 0;
			for (var i = 0; i < bytes.Length; i += 3)
			{
				if (maxB < bytes[i])
				{
					maxB = bytes[i];
				}
			}
			image.UnlockBits(bmpData);

			return maxB;
		}

		/// <summary>
		/// [molecule]
		/// 
		/// Calculates mean brightness for the given image.
		/// </summary>
		/// <param name="image">Input image.</param>
		/// <returns>Mean brightness value</returns>
		public static float CalculateMeanBrightness(Bitmap image)
		{
			if (image.PixelFormat == PixelFormat.Format24bppRgb
				|| image.PixelFormat == PixelFormat.Format32bppRgb
				|| image.PixelFormat == PixelFormat.Format32bppArgb)
			{
				return CalculateMeanBrightnessRGB(image);
			}

			return 0.0f;
		}

		/// <summary>
		/// [molecule]
		/// 
		/// Calculates mean brightness for the given color image with 24 or 32 bits per pixel.
		/// </summary>
		/// <param name="image">Input image.</param>
		/// <returns>Mean brightness value</returns>
		private static float CalculateMeanBrightnessRGB(Bitmap image)
		{
			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			BitmapData bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);

			int bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			uint sumR = 0, sumG = 0, sumB = 0;
			int width = bmpData.Width * 3,
				alignByte = bmpData.Stride - width;
			for (int i = 0; i < bytesSize;)
			{
				int limit = i + width;
				for (int j = i; j < limit; j += 3)
				{
					sumB += bytes[j];
					sumG += bytes[j + 1];
					sumR += bytes[j + 2];
				}
				i = limit + alignByte;
			}
			var res = sumR * Const.GRAYSCALE_R_COEFF + sumG * Const.GRAYSCALE_G_COEFF + sumB * Const.GRAYSCALE_B_COEFF;
			res /= image.Size.Width * image.Size.Height;

			image.UnlockBits(bmpData);

			return res;
		}

		/// <summary>
		/// [molecule]
		/// 
		/// Calculates frequencies of gray levels for the given image.
		/// Important: the image is supposed to be a 24-bit grayscale one!
		/// </summary>
		/// <param name="image">Input image.</param>
		/// <param name="freq">Array of frequencies.</param>
		public static void CalculateGrayFrequencies(Bitmap image, ref float[] freq)
		{
			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);

			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			freq = new float[Const.GRAYLEVELS_COUNT];
			for (var i = 0; i < bytesSize; i += 3)
			{
				freq[bytes[i]]++;
			}
			image.UnlockBits(bmpData);
		}

		public static int CalculateSumBrightness(Bitmap image)
		{
			var count = 0;

			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			for (var i = 0; i < bytes.Length; ++i)
			{
				count += bytes[i];
			}

			image.UnlockBits(bmpData);

			return count;
		}

		public static int CalculatePixels(Bitmap image, float val)
		{
			var count = 0;

			var rect = new Rectangle(0, 0, image.Size.Width, image.Size.Height);
			var bmpData = image.LockBits(rect, ImageLockMode.ReadWrite, PixelFormat.Format24bppRgb);
			var bytesSize = image.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy(bmpData.Scan0, bytes, 0, bytesSize);

			int width = bmpData.Width * 3,
				alignByte = bmpData.Stride - width;
			for (var i = 0; i < bytes.Length; )
			{
				var limit = i + width;
				for (var j = i; j < limit; j += 3)
				{
					var tempGray = bytes[j] * Const.GRAYSCALE_B_COEFF
									+ bytes[j + 1] * Const.GRAYSCALE_G_COEFF
									+ bytes[j + 2] * Const.GRAYSCALE_R_COEFF;
					if (tempGray > val)
					{
						++count;
					}
				}
				i = limit + alignByte;
			}

			image.UnlockBits(bmpData);

			return count;
		}

		public static float CalculateEntropy(float[] freq)
		{
			float sum = 0.0f;
			foreach (float val in freq) { sum += val; }

			float sum_1 = 1.0f/ sum;
			float res = 0;

			for (int i = 0; i < freq.Length; ++i)
			{
				if (freq[i] > 0)
				{
					float temp = freq[i] * sum_1;
					res -= (float)(temp * Math.Log(temp));
				}
			}

			return res;
		}
		#endregion
		#endregion

		#region - Harris Detector class. -
		/// <summary>
		/// Class to implement Harris detector.
		/// </summary>
		public class HarrisDetector
		{
			public const float DEFAULT_SIGMA = 2.0f;
			private const float TH_THRESHOLD = 1500.0f;

			/// <summary>
			/// Calculates markers denoting points found by the Harris detector for grayscales input image.
			/// </summary>
			/// <param name="image">Input image.</param>
			/// <param name="sigma">Harris-Laplace detector sigma value.</param>
			/// <returns>List of markers.</returns>
			public static List<CircleMarker> CalculateMarkers(Bitmap image, float sigma)
			{
				// calculate detector.
				var map = Calculate(image, sigma);
				return Generic.GetMarkers(map, 0.0f, (int)(sigma * 3));
			}

			/// <summary>
			/// Calculates markers denoting points found by the Harris detector for color input image.
			/// </summary>
			/// <param name="image">Input image.</param>
			/// <param name="sigma">Harris-Laplace detector sigma value.</param>
			/// <returns>List of markers.</returns>
			public static List<CircleMarker> CalculateColorMarkers(Bitmap image, float sigma)
			{
				// calculate detector.
				var map = CalculateColor(image, sigma);
				return Generic.GetMarkers(map, 0.0f, (int)(sigma*3));
			}

			/// <summary>
			/// Calculates 2D map for Harris detector responses for grayscale input image.
			/// </summary>
			/// <param name="image">Input image.</param>
			/// <param name="sigma">Harris-Laplace detector sigma value.</param>
			/// <returns>Resulting 2D map.</returns>
			public static float[,] Calculate(Bitmap image, float sigma)
			{
				// calculate  detector.
				var img = new Bitmap(image);
				ImageConversion.ConvertToGrayscale(img);

				//
				// calculate gradients.
				int width = img.Width, height = img.Height;
				var gradX = CalculateGradient(img, Orientation.Horizontal);
				var gradY = CalculateGradient(img, Orientation.Vertical);
				var gradXY = new float[height, width];

				// calculate [gradX^2], [gradY^2] and [gradX * gradY] matrices.
				for (int i = 0; i < height; ++i)
				{
					for (int j = 0; j < width; ++j)
					{
						gradXY[i, j] = gradX[i, j] * gradY[i, j];
						gradX[i, j] *= gradX[i, j];
						gradY[i, j] *= gradY[i, j];
					}
				}

				// perform smoothing of matrices.
				gradX = ImageFilter.CalculateGaussianMeans(gradX, sigma);
				gradY = ImageFilter.CalculateGaussianMeans(gradY, sigma);
				gradXY = ImageFilter.CalculateGaussianMeans(gradXY, sigma);

				var c = new float[height, width];
				const float k = 0.04f;
				for (int i = 0; i < height; ++i)
				{
					for (int j = 0; j < width; ++j)
					{
						float m11 = gradX[i, j];
						float m12 = gradXY[i, j];
						float m22 = gradY[i, j];

						float detM = m11 * m22 - m12 * m12;
						float traceM = m11 + m22;

						c[i, j] = detM - k * traceM * traceM;
						if (c[i, j] < TH_THRESHOLD) c[i, j] = 0.0f;
					}
				}
				c = ImageFilter.Max(c);

				return c;
			}

			/// <summary>
			/// Calculates 2D map for Harris detector responses for color input image.
			/// </summary>
			/// <param name="image">Input image.</param>
			/// <param name="sigma">Harris-Laplace detector sigma value.</param>
			/// <returns>Resulting 2D map.</returns>
			public static float[,] CalculateColor(Bitmap image, float sigma)
			{
				// calculate Harris detector.

				//
				// calculate gradients.
				const int colorChannels = 3;
				int width = image.Width, height = image.Height;
				var colors = new Bitmap[colorChannels];
				GetColorChannels(image, ref colors);

				var gradX = new List<float[,]>();
				var gradY = new List<float[,]>();
				var gradXY = new List<float[,]>();

				for (int i = 0; i < colorChannels; ++i)
				{
					gradX.Add(CalculateGradient(colors[i], Orientation.Horizontal));
					gradY.Add(CalculateGradient(colors[i], Orientation.Vertical));
					gradXY.Add(new float[height, width]);
				}

				// calculate [gradX^2], [gradY^2] and [gradX * gradY] matrices.
				var m11 = new float[height, width];
				var m22 = new float[height, width];
				var m12 = new float[height, width];
				for (int i = 0; i < height; ++i)
				{
					for (int j = 0; j < width; ++j)
					{
						gradXY[0][i,j] = gradX[0][i,j] * gradY[0][i,j];
						gradX[0][i,j] *= gradX[0][i,j];
						gradY[0][i,j] *= gradY[0][i,j];

						gradXY[1][i,j] = gradX[1][i,j] * gradY[1][i,j];
						gradX[1][i,j] *= gradX[1][i,j];
						gradY[1][i,j] *= gradY[1][i,j];

						gradXY[2][i,j] = gradX[2][i,j] * gradY[2][i,j];
						gradX[2][i,j] *= gradX[2][i,j];
						gradY[2][i,j] *= gradY[2][i,j];

						m11[i,j] = (gradX[0][i,j] + gradX[1][i,j] + gradX[2][i,j]);
						m22[i,j] = (gradY[0][i,j] + gradY[1][i,j] + gradY[2][i,j]);
						m12[i,j] = (gradXY[0][i,j] + gradXY[1][i,j] + gradXY[2][i,j]);
					}
				}

				// perform smoothing of matrices.
				m11 = ImageFilter.CalculateGaussianMeans(m11, sigma);
				m22 = ImageFilter.CalculateGaussianMeans(m22, sigma);
				m12 = ImageFilter.CalculateGaussianMeans(m12, sigma);

				var c = new float[height, width];
				const float k = 0.04f;
				for (int i = 0; i < height; ++i)
				{
					for (int j = 0; j < width; ++j)
					{
						float detM = m11[i, j] * m22[i, j] - m12[i, j] * m12[i, j];
						float traceM = m11[i, j] + m22[i, j];

						c[i,j] = detM - k * traceM * traceM;

						if (c[i, j] < TH_THRESHOLD) c[i, j] = 0.0f;
					}
				}
				c = ImageFilter.Max(c);

				return c;
			}

			/// <summary>
			/// Calculate points for Harris detector.
			/// </summary>
			/// <param name="image">Input image.</param>
			/// <param name="sigma">Harris-Laplace sigma value.</param>
			/// <returns>Set of corner points.</returns>
			public static List<Point> CalculateColorPoints(Bitmap image, float sigma)
			{
				// calculate Harris detector.
				var c = CalculateColor(image, sigma);
				return Generic.GetPoints(c, 0.0f);
			}
		}
		#endregion

		#region - Harris-Laplace Detector class. -
		/// <summary>
		/// Harris-Laplace detector as described in 
		/// 
		/// Mikolajczyk K., Schmid C. Indexing based on scale invariant interest points.
		/// 
		/// </summary>
		public class HarrisLaplaceDetector
		{
			private const int SCALES_COUNT = 6;
			private const float STARTING_SCALE = 1.5f;
			private const float SCALE_FACTOR = 1.4f;
			private const float TL_THRESHOLD = 10.0f;

			/// <summary>
			/// Finds corner points representing local maximums of the LoG response.
			/// </summary>
			/// <param name="harMaps">Set of Harris detector responses at different scale levels.</param>
			/// <param name="gaussMaps">Scale-space maps.</param>
			/// <param name="rads">Size of radiuses to define the marker size.</param>
			/// <returns>Set of circle markers storing location and radius of corner points.</returns>
			private static List<CircleMarker> ExtractSignificantPoints(List<List<Point>> harMaps, List<byte[,]> gaussMaps, float[] rads)
			{
				if (harMaps.Count == 0 || gaussMaps.Count == 0 || rads.Length == 0) { return null; }

				//
				// build additional map to refer to the level with maximal LoG response.
				var res = new List<CircleMarker>();

				// check all intermediate scales.
				for (int k = 1; k < SCALES_COUNT - 1; ++k)	// for each scale.
				{
					foreach (var point in harMaps[k])
					{
						int i = point.Y, j = point.X;
						if (gaussMaps[k][i, j] > gaussMaps[k - 1][i, j]
							&& gaussMaps[k][i, j] > gaussMaps[k + 1][i, j]
							&& gaussMaps[k][i, j] > TL_THRESHOLD)
						{
							res.Add(new CircleMarker(j, i, (int)rads[k]));
						}
					}
				}

				// check the smallest scale.
				foreach (var point in harMaps[0])
				{
					int i = point.Y, j = point.X;
					if (gaussMaps[0][i, j] > gaussMaps[1][i, j]
						&& gaussMaps[0][i, j] > TL_THRESHOLD)
					{
						res.Add(new CircleMarker(j, i, (int)rads[0]));
					}
				}

				// check the largest scale.
				const int idx = SCALES_COUNT - 1;
				const int idx_1 = idx - 1;
				foreach (var point in harMaps[idx])
				{
					int i = point.Y, j = point.X;
					if (gaussMaps[idx][i, j] > gaussMaps[idx_1][i, j]
						&& gaussMaps[idx][i, j] > TL_THRESHOLD)
					{
						res.Add(new CircleMarker(j, i, (int)rads[idx]));
					}
				}

				return res;
			}
			
			/// <summary>
			/// Calculate Harris-Laplace detector using iterative scaling to build-up a scale-space representation.
			/// This is a slightly faster method due to iterative building of scale-space maps. But it produces results
			/// which are different from that of the "standard" Harris-Laplace detector.
			/// </summary>
			/// <param name="img">Input image.</param>
			/// <returns>Set of circle markers storing location and radius of corner points.</returns>
			public static List<CircleMarker> CalculateFast(Bitmap img)
			{
				const float ksi = SCALE_FACTOR;	// scale change coefficient.
				const float scoef = 1.0f / SCALE_FACTOR;	// scale change coefficient for LoG operator.

				// Algorithm description.
				// 1. Calculate maps of different scales.
				// 2. Calculate Harris detector responses for each map.
				// 3. Remove points with low response values.
				// 4. Define scale for all remaining points.

				#region - Calculate corner points and LoG responses at different scales. -
				float curScale = STARTING_SCALE;
				var rads = new float[SCALES_COUNT];
				var harMaps = new List<List<Point>>();
				var gaussMaps = new List<byte[,]>();
				var scImg = new Bitmap(img);
				for (int i = 0; i < SCALES_COUNT; ++i, curScale *= ksi)
				{
					rads[i] = curScale * 3;

					scImg = ImageFilter.GaussianBlur(scImg, ksi);
					gaussMaps.Add(ImageConversion.ConvertToGrayscaleTable(scImg));

					// use Harris detector to select corner points at the current scale.
					var harMap = HarrisDetector.CalculateColorPoints(scImg, curScale * scoef);
					harMaps.Add(harMap);
				} 
				#endregion

				return ExtractSignificantPoints (harMaps, gaussMaps, rads);
			}

			/// <summary>
			/// Calculates Harris-Laplace detector as described in Mikolajczyk K., Schmid C. Indexing based on scale invariant interest points.
			/// </summary>
			/// <param name="img">Input image.</param>
			/// <returns>Set of circle markers storing location and radius of corner points.</returns>
			public static List<CircleMarker> Calculate(Bitmap img)
			{
				const float ksi = SCALE_FACTOR;	// scale change coefficient.
				const float scoef = 1.0f / SCALE_FACTOR;	// scale change coefficient for LoG operator.

				// Algorithm description.
				// 1. Calculate maps of different scales.
				// 2. Calculate Harris detector responses for each map.
				// 3. Remove points with low response values.
				// 4. Define scale for all remaining points.

				#region - Calculate corner points and LoG responses at different scales. -
				float curScale = STARTING_SCALE;
				var rads = new float[SCALES_COUNT];
				var harMaps = new List<List<Point>>();
				var gaussMaps = new List<byte[,]>();
				for (int i = 0; i < SCALES_COUNT; ++i, curScale *= ksi)
				{
					rads[i] = curScale * 3;

					// obtain image at the current scale.
					Bitmap scImg = ImageFilter.GaussianBlur(img, curScale);
					gaussMaps.Add(ImageConversion.ConvertToGrayscaleTable(scImg));

					// use Harris detector to select corner points at the current scale.
					var harMap = HarrisDetector.CalculateColorPoints(scImg, curScale * scoef);
					harMaps.Add(harMap);
				}
				#endregion

				return ExtractSignificantPoints(harMaps, gaussMaps, rads);
			}
		}
		#endregion

		#region - Image Evaluator class. -
		public class ImageEvaluator 
		{
			public static float EvaluateMR(Bitmap image)
			{
				var tempImg = new Bitmap(image);
				//var eval = 0.0f;

				var greyFreq = new float[Const.GRAYLEVELS_COUNT];

				// prepare auxiliary containers
				var width = image.Width;
				var height = image.Height;
				var square = width * height;
				var avLuma = CalculateMeanBrightness(image);

				// calculate grey frequencies
				ImageConversion.ConvertToGrayscale(tempImg);
				CalculateGrayFrequencies(tempImg, ref greyFreq);

				float _H = 0;
				for (var i = 0; i < greyFreq.Length; i++)
				{
					if (greyFreq[i] != 0)
					{
						greyFreq[i] /= square;
						_H -= (float)(greyFreq[i] * Math.Log(greyFreq[i]));// * log2_1;
					}
				}
				tempImg = (Bitmap)image.Clone();
				ImageFilter.EdgeDetector.Sobel(tempImg, (float)Math.Sqrt(avLuma));
				//tempImg.Save("edge.bmp");

				var _E = CalculateSumBrightness(tempImg) / 3.0f;
				var _Mu = CalculatePixels(tempImg, avLuma) / 3.0f;
				var eval = (float)(Math.Log(Math.Log(_E + Math.E)) * Math.Exp(_H) * _Mu / square);

				return eval;
			}

			public static float EvaluateMRMod(Bitmap image)
			{
				var tempImg = new Bitmap(image);
				var greyFreq = new float[Const.GRAYLEVELS_COUNT];

				// prepare auxiliary containers
				var width = image.Width;
				var height = image.Height;
				var square = width * height;
				var avLuma = CalculateMeanBrightness(image);

				// calculate grey frequencies
				ImageConversion.ConvertToGrayscale(tempImg);
				//tempImg.Save("gray.bmp");
				CalculateGrayFrequencies(tempImg, ref greyFreq);

				var _H = CalculateEntropy(greyFreq);

				tempImg = (Bitmap)image.Clone();
				ImageFilter.EdgeDetector.Sobel(tempImg, (float)Math.Sqrt(avLuma));
				//tempImg.Save("edge.bmp");

				var _E = CalculateSumBrightness(tempImg) / (3.0f * CalculatePixels(tempImg, 0));
				var _Mu = CalculatePixels(tempImg, avLuma) / 3.0f;
				var eval = (float)((256.0 - Math.Exp(_H)) / 192.0 + (square - _Mu) / (square) + (255.0 - _E) / 255.0);

				return eval;
			}
		}
		#endregion

		#region - Image sampling. -
		/// <summary>
		/// Predefined names for point detectors supported by [GetImagePoints] method.
		/// </summary>
		public static string[] POINT_DETECTORS = {
		                                           	"All Image Points",
		                                           	"Harris Detector",
		                                           	"Random Points (10%)",
		                                           	"Random Points (20%)",
		                                           	"Random Points (30%)",
		                                           	"Random Points (40%)",
		                                           	"Random Points (50%)",
		                                           	"Random Points (60%)",
		                                           	"Random Points (70%)",
		                                           	"Random Points (80%)",
		                                           	"Random Points (90%)"
		                                           };

		/// <summary>
		/// [molecule]
		/// 
		/// Returns selected points from the given image using specified point detector.
		/// If all image points are to be selected, then empty list is returned.
		/// </summary>
		/// <param name="img">Input image.</param>
		/// <param name="descr">Point detector name (see [POINT_DETECTORS] array).</param>
		/// <param name="rand">Random numbers generator.</param>
		/// <returns>List of selected points.</returns>
		public static List<Point> GetImagePoints (Bitmap img, string descr, Random rand)
		{
			#region - Use all points. -
			if (string.Compare(descr, POINT_DETECTORS[0]) == 0)
			{	// full description --> return empty set of points.
				return new List<Point>();
			} 
			#endregion

			#region - Use point detectors. -
			if (string.Compare(descr, POINT_DETECTORS[1]) == 0)
			{	// Harris detector.
				return HarrisDetector.CalculateColorPoints(img, HarrisDetector.DEFAULT_SIGMA);
			} 
			#endregion

			#region - Use random points. -
			int square = img.Height * img.Width;
			if (string.Compare(descr, POINT_DETECTORS[2]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.1 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[3]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.2 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[4]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.3 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[5]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.4 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[6]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.5 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[7]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.6 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[8]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.7 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[9]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.8 * square), rand);
			}
			if (string.Compare(descr, POINT_DETECTORS[10]) == 0)
			{	// Random sampling.
				return RandomImageSampler.CalculatePoints(img, (int)(0.9 * square), rand);
			}
			#endregion

			throw new Exception("[GetImagePoints]: Invalid argument exception. (" + descr + ")");
		}

		/// <summary>
		/// [molecule]
		/// 
		/// Returns subimage of the given image bounded by the specified rectangle.
		/// </summary>
		/// <param name="img">Source image.</param>
		/// <param name="argRect">Subimage bounding rectangle.</param>
		/// <returns>Requested subimage.</returns>
		public static Bitmap GetSubimage(Bitmap img, Rectangle argRect)
		{
			int width = argRect.Width, height = argRect.Height;
			//int stride = 3 * width + width % 4;

			//Rectangle rect = new Rectangle (0, 0, img.Size.Width, img.Size.Height);
			Rectangle rect = argRect;
			BitmapData bmpData = img.LockBits(rect, ImageLockMode.ReadOnly, PixelFormat.Format24bppRgb);
			int bytesSize = rect.Height * bmpData.Stride;
			var bytes = new byte[bytesSize];
			Marshal.Copy( bmpData.Scan0, bytes, 0, bytesSize);

			var newBytes = new byte[bytesSize];
			Array.Copy(bytes, newBytes, newBytes.Length);
			img.UnlockBits(bmpData);

			// copy new bytes into the resBmp.
			Rectangle rect2 = new Rectangle(0, 0, width, height);
			Bitmap resBmp = new Bitmap(width, height, PixelFormat.Format24bppRgb);
			BitmapData bmpData2 = resBmp.LockBits(rect2, ImageLockMode.WriteOnly, PixelFormat.Format24bppRgb);
			Marshal.Copy( newBytes, 0, bmpData2.Scan0, bytesSize);
			resBmp.UnlockBits(bmpData2);

			return resBmp;
		}

		#endregion

		#region - Random image sampling. -
		/// <summary>
		/// Class to sample points from the image.
		/// </summary>
		public class RandomImageSampler
		{
			/// <summary>
			/// [molecular]
			/// 
			/// Returns list of random pixels coordinates lying in the given image.
			/// </summary>
			/// <param name="img">Image to sample points from.</param>
			/// <param name="number">Number of points to sample.</param>
			/// <param name="rand">Random numbers generator.</param>
			/// <returns>Set of random points.</returns>
			public static List<Point> CalculatePoints(Bitmap img, int number, Random rand)
			{
				return StructMath.SampleRandomPoints(img.Size, number, rand);
			}
		}
		#endregion

		#region - EdgeOrientationMaskGenerator class. -
		/// <summary>
		/// [atomic]
		/// 
		/// Class to generate masks which are required to find lines with certain orientation.
		/// </summary>
		public class EdgeOrientationMaskGenerator
		{
			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateMasks8Lines (int size)
			{
				int size_2 = size/2;
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var mask0 = new float[size,size];
				for (int i = size_2; i < size; ++i ) {mask0[size_2, i] = 1.0f;}
				res.Add(mask0);

				//
				// create 45 degrees line mask.
				var mask45 = new float[size, size];
				for (int i = size_2; i < size; ++i)
				{
					int rowIdx = size - i - 1;
					mask45[rowIdx, i] = 1.0f;
				}
				res.Add(mask45);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				for (int i = 0; i < 3; ++i )
				{
					mask0 = MatrixMath.Rotate90CCW(mask0);
					mask45 = MatrixMath.Rotate90CCW(mask45);
					res.Add(mask0);
					res.Add(mask45);
				}

				return res;
			}

			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks with gaussian smoothing which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <param name="sigma">Gaussian smoothing parameter.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateMasks8Lines(int size, float sigma)
			{
				int size_2 = size / 2;
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var mask0 = new float[size, size];
				for (int i = size_2; i < size; ++i) { mask0[size_2, i] = 1.0f; }
				res.Add(mask0);

				//
				// create 45 degrees line mask.
				var mask45 = new float[size, size];
				for (int i = size_2; i < size; ++i)
				{
					int rowIdx = size - i - 1;
					mask45[rowIdx, i] = 1.0f;
				}
				res.Add(mask45);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				for (int i = 0; i < 3; ++i)
				{
					mask0 = MatrixMath.Rotate90CCW(mask0);
					mask45 = MatrixMath.Rotate90CCW(mask45);
					res.Add(mask0);
					res.Add(mask45);
				}

				//
				// Apply Gaussian smoothing.
				for (int i = 0; i < res.Count; ++i )
				{
					res[i] = ImageFilter.CalculateGaussianMeans(res[i], sigma);

					//
					// normalize mask so that sum of its elements = 1.
					var maxVal = MatrixMath.Sum(res[i]);
					res[i] = MatrixMath.Mul(res[i], 1.0f/maxVal);
				}

				return res;
			}

			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks using Gabor filter, which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateGaborMasks8LinesEx(int size)
			{
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var geParams = new GaussianEnvelopeParameters();
				geParams.Angle = 0;
				geParams.Center = new PointF(0.5f, 0);
				geParams.Magnitude = 1.0f;
				geParams.ScaleXSqr = 0.7f;
				geParams.ScaleYSqr = 6.0f;
				var mask0 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask0);
				var sum = MatrixMath.Sum(mask0);
				mask0 = MatrixMath.Mul(mask0, 1.0f / sum);
				res.Add(mask0);

				//
				// create 45 degrees line mask.
				geParams.Center = new PointF(0.5f, 0.5f);
				geParams.Angle = (float)(Math.PI / 4);
				var mask45 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask45);
				sum = MatrixMath.Sum(mask45);
				mask45 = MatrixMath.Mul(mask45, 1.0f / sum);
				res.Add(mask45);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				for (int i = 0; i < 3; ++i)
				{
					mask0 = MatrixMath.Rotate90CCW(mask0);
					mask45 = MatrixMath.Rotate90CCW(mask45);
					res.Add(mask0);
					res.Add(mask45);
				}

				return res;
			}

			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks using Gabor filter, which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateGaborMasks8Lines(int size)
			{
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var geParams = new GaussianEnvelopeParameters();
				geParams.Angle = 0;
				geParams.Center = new PointF(0.4f, 0);
				geParams.Magnitude = 1.0f;
				geParams.ScaleXSqr = 3.0f;
				geParams.ScaleYSqr = 3.0f;
				var mask0 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), (int)(size / 2), 1.5f, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask0);
				var sum = MatrixMath.Sum(mask0);
				mask0 = MatrixMath.Mul(mask0, 1.0f/sum);
				res.Add(mask0);

				//
				// create 45 degrees line mask.
				geParams.Center = new PointF(0.4f, 0.4f);
				//geParams.Angle = 0f;
				var mask45 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), (int)(size * 4 / 9), 1f, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask45);
				sum = MatrixMath.Sum(mask45);
				mask45 = MatrixMath.Mul(mask45, 1.0f / sum);
				res.Add(mask45);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				for (int i = 0; i < 3; ++i)
				{
					mask0 = MatrixMath.Rotate90CCW(mask0);
					mask45 = MatrixMath.Rotate90CCW(mask45);
					res.Add(mask0);
					res.Add(mask45);
				}

				return res;
			}
			
			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks using Gabor filter, which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateGaborMasks8LinesCentered(int size)
			{
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var geParams = new GaussianEnvelopeParameters();
				geParams.Angle = 0;
				geParams.Center = new PointF(0, 0);
				geParams.Magnitude = 1.0f;
				geParams.ScaleXSqr = 4.0f;
				geParams.ScaleYSqr = 4.0f;
				var mask0 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 2f, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask0);
				var sum = MatrixMath.Sum(mask0);
				mask0 = MatrixMath.Mul(mask0, 1.0f / sum);
				res.Add(mask0);

				//
				// create 22.5 degrees line mask.
				var mask225 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0.75f, 0.5f * (size - 5), 0, geParams);
				MatrixMath.NormalizeSelf(ref mask225);
				sum = MatrixMath.Sum(mask225);
				mask225 = MatrixMath.Mul(mask225, 1.0f / sum);
				res.Add(mask225);

				//
				// create 45 degrees line mask.
				var mask45 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 1, 0.5f*(size-3), 0, geParams);
				MatrixMath.NormalizeSelf(ref mask45);
				sum = MatrixMath.Sum(mask45);
				mask45 = MatrixMath.Mul(mask45, 1.0f / sum);
				res.Add(mask45);

				//
				// create 77.5 degrees line mask.
				var mask775 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 1, 0.5f * (size - 2), 0, geParams);
				MatrixMath.NormalizeSelf(ref mask775);
				sum = MatrixMath.Sum(mask775);
				mask775 = MatrixMath.Mul(mask775, 1.0f / sum);
				res.Add(mask775);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				var tempMask = MatrixMath.Rotate90CCW(mask0);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask225);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask45);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask775);
				res.Add(tempMask);

				return res;
			}

			/// <summary>
			/// [atomic]
			/// 
			/// Creates 8 square masks using Gabor filter, which represent 4 semi-axis lines and 4 semi-diagonal lines.
			/// </summary>
			/// <param name="size">Masks size.</param>
			/// <returns>Resulting list of masks.</returns>
			public static List<float[,]> CreateMasks8LinesCentered(int size)
			{
				var res = new List<float[,]>();

				//
				// create 0 degrees line mask.
				var geParams = new GaussianEnvelopeParameters();
				geParams.Angle = 0;
				geParams.Center = new PointF(0, 0);
				geParams.Magnitude = 1.0f;
				geParams.ScaleXSqr = 1.5f;
				geParams.ScaleYSqr = size;
				var mask0 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask0);
				var sum = MatrixMath.Sum(mask0);
				mask0 = MatrixMath.Mul(mask0, 1.0f / sum);
				res.Add(mask0);

				//
				// create 22.5 degrees line mask.
				geParams.Angle = (float)(Math.PI / 8);
				var mask225 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask225);
				sum = MatrixMath.Sum(mask225);
				mask225 = MatrixMath.Mul(mask225, 1.0f / sum);
				res.Add(mask225);

				//
				// create 45 degrees line mask.
				geParams.Angle = (float)(Math.PI / 4);
				var mask45 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask45);
				sum = MatrixMath.Sum(mask45);
				mask45 = MatrixMath.Mul(mask45, 1.0f / sum);
				res.Add(mask45);

				//
				// create 72.5 degrees line mask.
				geParams.Angle = (float)(Math.PI * 3f / 8);
				var mask775 = ImageFilter.GaborFilter.CreateFilterMaskScaleFree(new Size(size, size), 0, 0, 0, geParams);
				MatrixMath.NormalizeSelf(ref mask775);
				sum = MatrixMath.Sum(mask775);
				mask775 = MatrixMath.Mul(mask775, 1.0f / sum);
				res.Add(mask775);

				//
				// All other masks are obtained via rotation of mask0 and mask45.
				var tempMask = MatrixMath.Rotate90CCW(mask0);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask225);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask45);
				res.Add(tempMask);
				tempMask = MatrixMath.Rotate90CCW(mask775);
				res.Add(tempMask);

				return res;
			}
		}
		#endregion
	}
}
